package org.activequant.math.algorithms;

import org.apache.log4j.Logger;

/**
 * Fits sequence of 2-vectors {y_0, y_1, ... y_n}, {x_0, x_1, ... x_n} into a linear model:
 * <pre>
 * y = alpha + phi x
 * </pre>
 * where alpha is 2-vector, phi is 2x2 matrix.
 * <p>
 * Result is the learned matrix phi, vector alpha, and noise correlation matrix sigma.
 * <b>History:</b><br>
 *  - [10.02.2008] Created (Mike Kroutikov)<br>
 * Holds the following associated variables:
 * <ul>
 * <li>x0(EMAAccumulator)</li>
 * <li>x1(EMAAccumulator)</li>
 * <li>f0(EMAAccumulator)</li>
 * <li>f1(EMAAccumulator)</li>
 * <li>b00(EMAAccumulator)</li>
 * <li>b01(EMAAccumulator)</li>
 * <li>b10(EMAAccumulator)</li>
 * <li>b11(EMAAccumulator)</li>
 * <li>i00(EMAAccumulator)</li>
 * <li>i01(EMAAccumulator)</li>
 * <li>i10(EMAAccumulator)</li>
 * <li>i11(EMAAccumulator)</li>
 * <li>k00(EMAAccumulator)</li>
 * <li>k01(EMAAccumulator)</li>
 * <li>k10(EMAAccumulator)</li>
 * <li>k11(EMAAccumulator)</li>
 * <li>phi00(double)</li>
 * <li>phi01(double)</li>
 * <li>phi10(double)</li>
 * <li>phi11(double)</li>
 * <li>alpha0(double)</li>
 * <li>alpha1(double)</li>
 * <li>sigma00(double)</li>
 * <li>sigma01(double)</li>
 * <li>sigma10(double)</li>
 * <li>sigma11(double)</li>
 * <li>solved(boolean)</li>
 * <li>EPSILON(double)</li>
 * <li>muSolved(boolean)</li>
 * <li>mu0(double)</li>
 * <li>mu1(double)</li>
 * </ul>
 *  @author Mike Kroutikov
 */
public class EMARegression2 {
	
	private final Logger log = Logger.getLogger(getClass());
	/**
	 * private final EMAAccumulator x0 = new EMAAccumulator();
	 */
	private final EMAAccumulator x0 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator x1 = new EMAAccumulator();
	 */
	private final EMAAccumulator x1 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator f0 = new EMAAccumulator();
	 */
	private final EMAAccumulator f0 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator f1 = new EMAAccumulator();
	 */
	private final EMAAccumulator f1 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator b00 = new EMAAccumulator();
	 */
	private final EMAAccumulator b00 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator b01 = new EMAAccumulator();
	 */
	private final EMAAccumulator b01 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator b10 = new EMAAccumulator();
	 */
	private final EMAAccumulator b10 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator b11 = new EMAAccumulator();
	 */
	private final EMAAccumulator b11 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator i00 = new EMAAccumulator();
	 */
	private final EMAAccumulator i00 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator i01 = new EMAAccumulator();
	 */
	private final EMAAccumulator i01 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator i10 = new EMAAccumulator();
	 */
	private final EMAAccumulator i10 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator i11 = new EMAAccumulator();
	 */
	private final EMAAccumulator i11 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator k00 = new EMAAccumulator();
	 */
	private final EMAAccumulator k00 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator k01 = new EMAAccumulator();
	 */
	private final EMAAccumulator k01 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator k10 = new EMAAccumulator();
	 */
	private final EMAAccumulator k10 = new EMAAccumulator();
	/**
	 * private final EMAAccumulator k11 = new EMAAccumulator();
	 */
	private final EMAAccumulator k11 = new EMAAccumulator();
	/**
	 * Constructs an EMARegression2 using the given period(int) to set the periods and lambdas of its associated EMAAccumulators
	 * @param period
	 */
	public EMARegression2(int period) {
		x0.setPeriod(period);
		x1.setPeriod(period);
		f0.setPeriod(period);
		f1.setPeriod(period);
		b00.setPeriod(period);
		b01.setPeriod(period);
		b10.setPeriod(period);
		b11.setPeriod(period);
		i00.setPeriod(period);
		i01.setPeriod(period);
		i10.setPeriod(period);
		i11.setPeriod(period);
		k00.setPeriod(period);
		k01.setPeriod(period);
		k10.setPeriod(period);
		k11.setPeriod(period);
	}
	/**
	 * Accumulates the given yy0(double), yy1(double), xx0(double) and xx1(double)
	 * @param yy0
	 * @param yy1
	 * @param xx0
	 * @param xx1
	 */
	public void accumulate(double yy0, double yy1, double xx0, double xx1) {
		double ff0 = yy0 - xx0;
		double ff1 = yy1 - xx1;
		
		double x0prev = x0.getMeanValue();
		x0.accumulate(xx0);
		double x1prev = x1.getMeanValue();
		x1.accumulate(xx1);

		double f0prev = f0.getMeanValue();
		f0.accumulate(ff0);
		double f1prev = f1.getMeanValue();
		f1.accumulate(ff1);
		
		if(log.isDebugEnabled()) {
			log.debug("xx0=" + (xx0 - x0prev));
			log.debug("xx1=" + (xx1 - x1prev));
			log.debug("ff0=" + (ff0 - f0prev));
			log.debug("ff1=" + (ff1 - f1prev));
		}

		b00.accumulate((xx0 - x0prev) * (xx0 - x0.getMeanValue()));
		b01.accumulate((xx0 - x0prev) * (xx1 - x1.getMeanValue()));
		b10.accumulate((xx1 - x1prev) * (xx0 - x0.getMeanValue()));
		b11.accumulate((xx1 - x1prev) * (xx1 - x1.getMeanValue()));

		i00.accumulate((ff0 - f0prev) * (ff0 - f0.getMeanValue()));
		i01.accumulate((ff0 - f0prev) * (ff1 - f1.getMeanValue()));
		i10.accumulate((ff1 - f1prev) * (ff0 - f0.getMeanValue()));
		i11.accumulate((ff1 - f1prev) * (ff1 - f1.getMeanValue()));

		k00.accumulate((ff0 - f0prev) * (xx0 - x0.getMeanValue()));
		k01.accumulate((ff0 - f0prev) * (xx1 - x1.getMeanValue()));
		k10.accumulate((ff1 - f1prev) * (xx0 - x0.getMeanValue()));
		k11.accumulate((ff1 - f1prev) * (xx1 - x1.getMeanValue()));
		
		if(log.isDebugEnabled()) {
			log.debug("x=" + x0.getMeanValue() + " " + x1.getMeanValue());
			log.debug("f=" + f0.getMeanValue() + " " + f1.getMeanValue());
			log.debug("b=" + b00.getMeanValue() + " " + b01.getMeanValue() + " " + b10.getMeanValue() + " " + b11.getMeanValue());
			log.debug("k=" + k00.getMeanValue() + " " + k01.getMeanValue() + " " + k10.getMeanValue() + " " + k11.getMeanValue());
		}

		solved = false;
	}
	/**
	 * private double phi00 = 0.0;
	 */
	private double phi00 = 0.0;
	/**
	 * private double phi01 = 0.0;
	 */
	private double phi01 = 0.0;
	/**
	 * private double phi10 = 0.0;
	 */
	private double phi10 = 0.0;
	/**
	 * private double phi11 = 0.0;
	 */
	private double phi11 = 0.0;
	/**
	 * private double alpha0 = 0.0;
	 */
	private double alpha0 = 0.0;
	/**
	 * private double alpha1 = 0.0;
	 */
	private double alpha1 = 0.0;
	/**
	 * private double sigma00 = 0.0;
	 */
	private double sigma00 = 0.0;
	/**
	 * private double sigma01 = 0.0;
	 */
	private double sigma01 = 0.0;
	/**
	 * private double sigma10 = 0.0;
	 */
	private double sigma10 = 0.0;
	/**
	 * private double sigma11 = 0.0;
	 */
	private double sigma11 = 0.0;
	/**
	 * private boolean solved = false;
	 */
	private boolean solved = false;
	/**
	 * private static final double EPSILON = 1.e-10;
	 */
	private static final double EPSILON = 1.e-10;
	/**
	 * Returns whether <strong>detB > EPSILON*detK && detK > EPSILON*detB</strong>, where <br/>
	 * <code>detB = b00.meanValue * b11.meanValue -  b01.meanValue * b10.meanValue</code> and<br/>
	 * <code>detK = abs(k00.meanValue * k11.meanValue -  k01.meanValue * k10.meanValue)</code>
	 * @return
	 */
	public boolean isReady() {
		double detB = b00.getMeanValue() * b11.getMeanValue() - b01.getMeanValue() * b10.getMeanValue();
		double detK = k00.getMeanValue() * k11.getMeanValue() - k01.getMeanValue() * k10.getMeanValue();
		detK = Math.abs(detK);
		
		return detB > EPSILON * detK && detK > EPSILON * detB;
	}
	/**
	 * Sets the associated solved(boolean) to true, and the associated muSolved(boolean) to false.<br/>
	 * recalculates the associated phi00(double), phi01(double), phi10(double), phi11(double), alpha0(double), alpha1(double), sigma00(double), sigma01(double), sigma10(double), sigma11(double)
	 */
	private void solve() {
		solved = true;
		muSolved = false;
		
		double b00 = this.b00.getMeanValue();
		double b01 = this.b01.getMeanValue();
		double b10 = this.b10.getMeanValue();
		double b11 = this.b11.getMeanValue();
		
		double detB = b00 * b11 - b01 * b10;
		double detK = k00.getMeanValue() * k11.getMeanValue() - k01.getMeanValue() * k10.getMeanValue();
		
		if(detB < EPSILON * Math.abs(detK)) return;
		
		// reverse the matrix B
		double revB00 =   b11 / detB;
		double revB01 = - b10 / detB;
		double revB10 = - b01 / detB;
		double revB11 =   b00 / detB;
		
		phi00 = k00.getMeanValue() * revB00 + k01.getMeanValue() * revB10;
		phi01 = k00.getMeanValue() * revB01 + k01.getMeanValue() * revB11;
		phi10 = k10.getMeanValue() * revB00 + k11.getMeanValue() * revB10;
		phi11 = k10.getMeanValue() * revB01 + k11.getMeanValue() * revB11;
		
		alpha0 = f0.getMeanValue() - phi00 * x0.getMeanValue() - phi01 * x1.getMeanValue();
		alpha1 = f1.getMeanValue() - phi10 * x0.getMeanValue() - phi11 * x1.getMeanValue();

		
		double n00 = k00.getMeanValue() * revB00 + k01.getMeanValue() * revB10;
		double n01 = k00.getMeanValue() * revB01 + k01.getMeanValue() * revB11;
		double n10 = k10.getMeanValue() * revB00 + k11.getMeanValue() * revB10;
		double n11 = k10.getMeanValue() * revB01 + k11.getMeanValue() * revB11;
		
		sigma00 = i00.getMeanValue();
		sigma01 = i01.getMeanValue();
		sigma10 = i10.getMeanValue();
		sigma11 = i11.getMeanValue();
		
		sigma00 -= n00 * k00.getMeanValue() + n01 * k01.getMeanValue();
		sigma01 -= n00 * k10.getMeanValue() + n01 * k11.getMeanValue();
		sigma10 -= n10 * k00.getMeanValue() + n11 * k01.getMeanValue();
		sigma11 -= n10 * k10.getMeanValue() + n11 * k11.getMeanValue();
		
//		log.info("phi=[" + (phi00 + 1.) + " " + phi01 + " " + phi10 + " " + (phi11 + 1.));
//		log.info("alpha=[" + alpha0 + " " + alpha1);
//		log.info("sigma=[" + sigma00 + " " + sigma01 + " " + sigma10 + " " + sigma11);
	}
	/**
	 * if the associated solve(boolean) is false, it invokes <code>solve()</code>.<br/>
	 * returns a double[] array holding <strong>{phi00+1, phi01, phi10, phi11+1}</strong>
	 * @return
	 */
	public double [] getPhi() {
		if(!solved) solve();
		return new double [] {phi00 + 1., phi01, phi10, phi11 + 1.};
	}
	/**
	 * if the associated solve(boolean) is false, it invokes <code>solve()</code>.<br/>
	 * returns a double[] array holding <strong>{alpha0, alpha1}</strong>
	 * @return
	 */
	public double [] getAlpha() {
		if(!solved) solve();
		return new double [] {alpha0, alpha1};
	}
	/**
	 * if the associated solve(boolean) is false, it invokes <code>solve()</code>.<br/>
	 * if the associated muSolved(boolean) is false, it invokes <code>muSolved()</code>.<br/>
	 * returns a double[] array holding <strong>{mu0, mu1}</strong>
	 * @return
	 */
	public double [] getMu() {
		if(!solved) solve();
		if(!muSolved) muSolve();
		
		return new double [] {mu0, mu1};
	}
	/**
	 * private boolean muSolved = false;
	 */
	private boolean muSolved = false;
	/**
	 * private double mu0 = 0.0;
	 */
	private double mu0 = 0.0;
	/**
	 * private double mu1 = 0.0;
	 */
	private double mu1 = 0.0;
	/**
	 * recalculates the associated mu0(double) and mu1(double), and sets the associated muSolved(boolean) to true
	 */
	private void muSolve() {
		double detB = b00.getMeanValue() * b11.getMeanValue() - b01.getMeanValue() * b10.getMeanValue();
		double detK = k00.getMeanValue() * k11.getMeanValue() - k01.getMeanValue() * k10.getMeanValue();
		
		if(Math.abs(detK) < EPSILON * detB) return;

		double revK00 =   k11.getMeanValue() / detK;
		double revK01 = - k10.getMeanValue() / detK;
		double revK10 = - k01.getMeanValue() / detK;
		double revK11 =   k00.getMeanValue() / detK;
		
		double bk00 = b00.getMeanValue() * revK00 + b01.getMeanValue() * revK10;
		double bk01 = b00.getMeanValue() * revK01 + b01.getMeanValue() * revK11;
		double bk10 = b10.getMeanValue() * revK00 + b11.getMeanValue() * revK10;
		double bk11 = b10.getMeanValue() * revK01 + b11.getMeanValue() * revK11;
	
		mu0 = x0.getMeanValue() - bk00 * f0.getMeanValue() - bk01 * f1.getMeanValue();
		mu1 = x1.getMeanValue() - bk10 * f0.getMeanValue() - bk11 * f1.getMeanValue();
		
		muSolved = true;
	}
	/**
	 * if the associated solve(boolean) is false, it invokes <code>solve()</code>.<br/>
	 * returns a double[] array holding <strong>{sigma00, sigma01, sigma10, sigma11}</strong>
	 * @return
	 */
	public double [] getSigma() {
		if(!solved) solve();
		return new double [] {sigma00, sigma01, sigma10, sigma11};
	}

	public String toString() {
		if(!solved) solve();
		
		return "EMARegresion2: period=" + x0.getPeriod() 
			+ ", numSamples=" + x0.getNumSamples()
			+ ", phi=[" + (phi00 + 1.) + " " + phi01 + " " + phi10 + " " + (phi11 + 1.) + "]" 
			+ ", alpha=[" + alpha0 + " " + alpha1 + "]" 
			+ ", sigma=[" + sigma00 + " " + sigma01 + " " + sigma10 + " " + sigma11 + "]";
	}
}
